%Written by Michael R. Springborn
%This version: Dec. 16, 2011

% Free to use and modify with proper acknowledgement. No warranties. 
%*******************************************************************************
%                         ***GENERAL DESCRIPTION***     

%This code analyzes and summarizes output from numerical simulation of the 
%Bayesian posterior (generated runBayesEst.m).  Results here are reported in Tables
%3 and 4 of Lieli and Springborn for "Bayes (logit)" and "Bayes (cauchit)"
%*******************************************************************************

clear; clc
work_dir='D:\Bayesian_Analysis'
cd(work_dir)

%% Load training data set
importfile('D:\data\data_WRA.csv'); %For descriptions of field names see Excel file of same name in same folder.

%Select the link model for modeling binomial outcomes: 
model='logit';
%model='cauchit';

%Assumed base rate of weeds (unconditional probability that Y=1).  Note: referred to in manuscript as "tau".
Q=0.05; 
%Q=0.02 

%Economic cutoff params
ben=141000; %benefit
dam=2068000; %damage
r=0.03; %discount rate
c_const=0.0682; %constant cutoff

%custom cutoff using undesireable traits
dam_c = inline('928300+786800.*Sc_Undes');
c_udt=ben./(ben+dam_c(Sc_Undes)-ben);

%MCMC parameters
burn=1000; %burn-in period to discard from each chain
s=20000; %iterations in each chain
M=11; %Number of chains

%SIMP: Model for the density of the covariate sample, f(x;lambda) in Equation (13) of Lieli and Springborn.  [Also see Appendix B in Lieli and Springborn]
   %Reported results in Lieli and Springborn are for simp=1 (mean and variance of x set to their maximum likelihood estimates)
   %See "CosGibbsMH.m" for detail on all four variations of increasing complexity.
SIMP=1:4; 

%Initialize arrays for collecting statistics:
Unb1=cell(length(Q),max(SIMP),M);
Unb0=Unb1; Pdist=Unb1; Ep=Unb1; ban=Unb1; ban_c=Unb1; 
sens=nan*ones(length(Q),max(SIMP),M); spec=sens; sens_c=sens; spec_c=sens; Eben_c=sens; Eben=sens;Eben_opendoor=sens;

%Prefix of location where MCMC output is stored
dest=['D:\results\mcmc_univariate\Q']; 

for d=1:length(Q)
    Qd=Q(d);
    for simp=1:4 %different models for dealing with distribution of covariate
        for m=1:M %number of chains
            S=load([dest num2str(Qd) '\simp' num2str(simp) '\' model '\Q' num2str(Qd) '_simp' num2str(simp) '_m' num2str(m) '_s' num2str(s) '.mat']);
            X=S.X; y=S.y; ybar=S.ybar; N=size(y,1);
            unind=find(S.beta1s~=[0;S.beta1s(1:end-1)]); %find index of where chain jumps to a new value
            Unb1{d,simp,m}=S.beta1s(unind(unind>burn)); %...exclude draws from burn in period and duplicated draws.
            Unb0{d,simp,m}=S.beta0s(unind(unind>burn));  
            %calculate fitted probabilities
            XB=[ones(size(X,1),1) X]*[Unb0{d,simp,m}' ; Unb1{d,simp,m}'];
            if strcmp(model,'logit')
                Pdist{d,simp,m}=1./(1+exp(-XB));
            elseif strcmp(model,'cauchit')
                Pdist{d,simp,m}=.5 + atan(XB)/pi;
            end
            %calculated desired stats
                %No suffix --> constant cutoff case
                %_c suffix --> custom cutoff case
            Ep{d,simp,m}   =mean(Pdist{d,simp,m},2);
            ban{d,simp,m}  =Ep{d,simp,m}>c_const; %Bayes bans
            ban_c{d,simp,m}=Ep{d,simp,m}>c_udt; 
            sens(d,simp,m) = sum(ban{d,simp,m}(y==1))/sum(y); %Sensitivity 
            spec(d,simp,m) = sum(1-ban{d,simp,m}(y==0))/sum(1-y); %Specificity
            Eben(d,simp,m) = Qd*(1-sens(d,simp,m))*(ben-dam)/r + (1-Qd)*spec(d,simp,m)*ben/r; %Expected benefits
            Eben_opendoor(d,simp,m) = -Qd*sens(d,simp,m)*(ben-dam)/r - (1-Qd)*(1-spec(d,simp,m))*ben/r; %If we have open door baseline...
            sens_c(d,simp,m) = sum(ban_c{d,simp,m}(y==1))/sum(y); % 
            spec_c(d,simp,m) = sum(1-ban_c{d,simp,m}(y==0))/sum(1-y); %
                accept_c=find(ban_c{d,simp,m}==0); %accepteds
                acceptweed_c= ban_c{d,simp,m}==0 & y==1;
                acptwdloss=((ben-dam_c(Sc_Undes))/r)'*acceptweed_c;
                acptnwdgain=ben/r*sum(1-y(accept_c));
            Eben_c(d,simp,m)=((Qd/ybar)*acptwdloss + ((1-Qd)/(1-ybar))*acptnwdgain)/N; %
        end
        display(['Simp = ' num2str(simp)])
    end
end

%% Display results:
%Dimensions for stats are {Qd,simp,m}. 

format short
%For table 3 in Lieli and Springborn
%d=1;
d=2;
squeeze(Eben(d,:,:))
squeeze(sens(d,:,:))
squeeze(spec(d,:,:))

squeeze(Eben_c(d,:,:)) %For table 4 in Lieli and Springborn
squeeze(sens_c(d,:,:))
squeeze(spec_c(d,:,:))

%squeeze(Eben_opendoor(1,:,:)) 
%squeeze(Eben_opendoor(2,:,:))

return

